
library(tidyverse)
library(rio)
library(cregg)
library(patchwork)

# set replication folder as working directory
setwd("~replication")

load("data_genderedcost_long.rdata")

# only include completed answers - and answers given before deadline
# 2021-12-20 21:49:52 was the last response within the time frame
df_long <- df_long %>% 
  filter(SurveyStatus==2)

df_long <- df_long %>% 
  filter(SurveyEndTime<="2021-12-20 21:49:52")

##### FIGURE 5 
## SUBSET WOMEN AND MEN's PREFERENCES FOR INTERAKTION BETWEEN POWER AND QUALITY OF WORKING ENVIRONMENT
# create interaction as column
df_long$inter_position_environment <-
  interaction(df_long$position, df_long$work_environment, sep = " and ")

# subset on men and women
df_men <- df_long %>% 
  dplyr::filter(sex=="Man")

mm_men_pos_env <- cregg::mm(df_men, # data
                            choice~inter_position_environment, # formula
                            id = id)

mm_men_pos_env$sex <- "Men" # add sex variable

df_women <- df_long %>% 
  dplyr::filter(sex=="Woman")

mm_women_pos_env <- cregg::mm(df_women, # data
                              choice~inter_position_environment, # formula
                              id = id)

mm_women_pos_env$sex <- "Women" # add sex variable

#combine estimates from the two subsets, add attribute variable for color and rearrange order of levels
mm_sex_subsets <- bind_rows(mm_men_pos_env, mm_women_pos_env) %>% 
  mutate(sex = factor(sex)) %>%
  mutate(position = case_when(str_detect(level, "Common")~"Common member",
                              str_detect(level, "large")~"Chair, large influence",
                              str_detect(level, "small")~"Chair, small influence"),
         environment = case_when(str_detect(level, "respect")~"Characterized by mutual respect",
                                 str_detect(level, "sexism")~"Several experienced \nsexism and harassment",
                                 TRUE~"Several experienced harassment")) %>%
  mutate(environment = factor(environment, levels = c("Characterized by mutual respect",
                                                      "Several experienced harassment",
                                                      "Several experienced \nsexism and harassment")))


#### SUBSET MM FOR INTERACTION BETWEEN POSITION AND QUALITY
## LEFT HAND-SIDE OF FIGURE 5
mm_pos_env_plot <- mm_sex_subsets %>% 
  ggplot(data=., aes(y = position, x = estimate, color = sex, shape = sex)) +
  geom_point(aes(shape = sex, color = sex), position = position_dodge2(width = 0.5)) +
  geom_linerange(aes(xmin=estimate-std.error*1.39,
                     xmax=estimate+std.error*1.39),
                 position = position_dodge2(width = 0.5)) +
  xlab("Marginal mean") +
  ylab("") +
  scale_color_grey("") +
  scale_shape_discrete("") +
  scale_x_continuous(breaks = seq(0.2,0.9,0.05), labels = seq(0.2,0.9,0.05)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  theme_bw() +
  facet_wrap(~environment, scales = "free_y", ncol = 1) +
  theme(legend.position = "top", panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())


######### GROUP DIFFERENCES IN MARGINAL MEANS
### Interaction: Three-way interaction: position x environment quality x gender
mm_int_pos_env <- mm_diffs(df_long, # data
                           choice~inter_position_environment, # formula
                           ~sex)


mm_int_pos_env <- mm_int_pos_env %>% 
  mutate(position = case_when(str_detect(level, "Common")~"Common member",
                              str_detect(level, "large")~"Chair, large influence",
                              str_detect(level, "small")~"Chair, small influence"),
         environment = case_when(str_detect(level, "respect")~"Characterized by \nmutual respect",
                                 str_detect(level, "sexism")~"Several experienced \nsexism and harassment",
                                 TRUE~"Several experienced \nharassment")) %>% 
  mutate(level = case_when(str_detect(level, "respect")~"Characterized by \nmutual respect",
                           str_detect(level, "sexism")~"Several experienced \nsexism and harassment",
                           TRUE~"Several experienced \nharassment"))

#### RIGHT HAND SIDE OF FIGURE 5
mm_int_pos_env_plot <- mm_int_pos_env %>% 
  ggplot(data=., aes(y = position, x = estimate)) +
  geom_point(position = position_dodge2(width = 0.25)) +
  geom_linerange(aes(xmin=lower,
                     xmax=upper),
                 position = position_dodge2(width = 0.25)) +
  xlab("Differences\n(women - men)") +
  ylab("") +
  scale_color_grey("") +
  scale_shape_discrete("") +
  scale_x_continuous(labels = seq(-0.2,0.2,0.05), breaks = seq(-0.2,0.2,0.05)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  facet_wrap(~environment, ncol = 1) +
  theme(legend.position = "none", panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())


### one big plot
mm_pos_env_plot + mm_int_pos_env_plot + plot_layout(widths = c(3, 1))

ggsave("figure5.pdf", width = 8, height = 4)
ggsave("figure5.png", width = 8, height = 4, dpi = 600)
